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We present an analysis of the work performed on a system of interest that is kept thermally 
isolated during the switching of a control parameter. We show that there exists, for a certain class 
of systems, a finite-time family of switching protocols for which the work is equal to the quasistatic 
value. These optimal paths are obtained within linear response for systems initially prepared in a 
canonical distribution. According to our approach, such protocols are composed of a linear part 
plus a function which is odd with respect to time reversal. For systems with one degree of freedom, 
we claim that these optimal paths may also lead to the conservation of the corresponding adiabatic 
invariant. This points to an interesting connection between work and the conservation of the volume 
enclosed by the energy shell. To illustrate our findings, we solve analytically the harmonic oscillator 
and present numerical results for certain anharmonic examples. 
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I. INTRODUCTION 

The accumulated and organized knowledge that we call 
thermodynamics has been one of the main pillars of our 
physical understanding of the world around us. However, 
the only processes that are fully describable by means of 
classical thermodynamics are quasistatic ones, i.e., pro¬ 
cesses which are a succession of equilibrium states [1]. 
On the other hand, real thermodynamic processes hap¬ 
pen in finite time [2-4] and hence drive the system out of 
equilibrium. In this case, the second law imposes certain 
limits to the energy exchange between a system of inter¬ 
est and an external agent. How to get as close as possible 
to the minimal energetic cost of driving a system from 
one state to another in finite time remains then a crucial 
question. Thus, it is very desirable to develop a general 
method to solve such optimization problem. 

Thermodynamic processes can be performed under dif¬ 
ferent constraints. The system of interest can be kept, for 
instance, in contact with a heat bath during the time in¬ 
terval its externally controlled parameter is switched. In 
this situation, the minimal energetic cost is equal to the 
difference of Helmholtz free energies. Thereby, one of the 
many applications of such optimal finite-time processes is 
the estimation of free-energy differences [5-10]. A major 
breakthrough in this problem was achieved by Jarzynski 
[11] and Crooks [12]. Through their results, finite-time 
processes can be used without leading to a biased esti¬ 
mation [13, 14]. Nevertheless, one needs to sample ex¬ 
tremely rare events in order to have reliable estimates 
[15]. Besides, different demands for better efficiency in 
finite time have increased the interest on optimal control 
of thermodynamic systems [16-24]. 

There are, at the moment, two main ways of finding op¬ 
timal finite-time processes under isothermal conditions: 
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stochastic models [25-27] and linear response theory [28- 
31]. The results obtained so far within the stochastic ap¬ 
proach show intriguing, interesting, and not well under¬ 
stood features which appear only for sufficiently fast pro¬ 
cesses. On the other hand, the linear response approach 
provides an analytical treatment of a broader class of sys¬ 
tems, although limited to quasiequilibrium processes. In 
the context of thermally isolated systems, the analysis of 
optimal paths has followed along the same lines as the 
stochastic approach for isothermal processes [32]. How¬ 
ever, an analytical description through stochastic meth¬ 
ods is very restricted to linear systems and leaves open 
questions about what happens in the nonlinear case. 

In the present work, we study the problem of finding 
optimal finite-time processes in thermally isolated sys¬ 
tems via linear response theory. We focus on the regime 
in which the variation of the externally controlled pa¬ 
rameter is small but has arbitrary speed. In Sec. H, we 
derive an expression for the excess work [29, 33, 34], i.e., 
a quantity that characterizes the energetic cost along a 
given process. In Secs. HI and IV, this expression is em¬ 
ployed to explain the numerical results of simple linear 
and nonlinear systems. The results obtained show unex¬ 
pected features from the point of view of usual thermo¬ 
dynamic wisdom. In Sec. V, we connect the excess work 
to the adiabatic invariant, suggesting that every time the 
former vanishes the latter is conserved. We summarize 
and conclude in Secs. VI and VH. 


II. THERMALLY ISOLATED SYSTEMS AND 
THE EXCESS WORK 

Let us consider the following setup: first, we keep our 
system of interest in contact with a heat bath until its 
relaxation to the Boltzmann-Gibbs distribution, 

Pe,(r; Ao) = exp(-/37^(r; Ao))/Z(/?, Aq) , (1) 

where F is a point in phase space, /3 = (fc^T)”^, with T 
and ks being the temperature of the heat bath and the 
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FIG. 1. (Color online) Schematic representation of a given 
protocol X{t) performed by the external agent while the sys¬ 
tem is kept thermally isolated. 

Boltzmann constant, respectively. The quantity Z{/3, Aq) 
denotes the partition function of the system given by 
Z{(3,Xq) = f dr exp(—/3H(r; Ao)) and Aq is the initial 
value of our control parameter A. Second, the system is 
decoupled from the reservoir and kept thermally isolated 
while the external agent switches A from Aq to Xf ac¬ 
cording to a given protocol (see Fig. 1). We express the 
protocol X(t) as follows: 

X(t) = Ao -f i5A g(t ), (2) 


Therefore, the nonequilibrium average of the generalized 
force can be calculated by means of linear response theory 
[36, 37]. It reads 

= +X^SXg{t) - SX dsMt- s)g{s), 

( 6 ) 

where (•)o denotes an average on the initial ensemble, 
given by Eq. (1), and the subscript refers to the value 
Aq. The second term in the right-hand side of Eq. (6) 
describes the instantaneous response, which is due to 
8%/dX being a function of the external control A [36, 37] . 
In particular, we have 

The second term describes the delayed response and (po (t) 
is the so-called response function. It will be convenient 
to express it in terms of the relaxation function, d'o(t) 
[36, 37]. This can be done as follows: 

Mt) = -^it) = -l3j^{Co{t)-C), ( 8 ) 

where C'o(t) = (A(0)Al(t))o is the correlation function of 
A = d'H/dX and the constant C is given by [37] 

pOO 

C = lim € dt e~'^^ Co{t). (9) 

<^-*■0 Jo 


where g{t) is such that g{to) = 0 and g{tf) = 1. Thus, 
the variation in A in the time interval r = tf — to is 
()A = Xf — Ao. 

If we consider just a single realization of the protocol 
X{t), the probability of reaching the value W of the work 
performed is given by V{W) = (5(>V — W[r(])), where 
Ft is a given trajectory in phase space [35]. This means 
that ViyV) can be calculated from an average over all the 
possible trajectories or realizations [11]. The thermody¬ 
namic work then reads 

pCO 

w= dWV{yV)W, (3) 


or, equivalently. 


W = 



dt 


dt dX 


(4) 


where A denotes the nonequilibrium average of the ob¬ 
servable A. 

Assuming that j(5A5(t)/Aol ^ 1 for Iq < ^ < tf, we can 
treat the effects of the generalized force d'H/dX pertur- 
batively. In other words, once our system of interest is 
described by a Hamiltonian 'H[X{t)], we can expand it in 
powers of (5A as follows: 


n[xit)] = n{Xo) + sxg{t) 


m 

dX 


0 ( 2 ). 


(5) 


Therefore, Eq. (6) can be rewritten after an integration 
by parts as 


fw = (f 


^A 


/■ 


*0 

»„(„)- 


( 10 ) 


where Ifo = 'I'o(O) — Einally, substituting Eq. (10) 
in expression (4), we obtain 


+ (a)2 I"' dt ^ dt' (11) 

using the boundary conditions for g{t). The first two 
terms of the previous expression do not depend on the 
protocol g{t). Indeed, it can be verihed (see Appendix 
B) that they are the first terms of the series expansion of 
the quasistatic work for ^A/Aq 1. The last term clearly 
depends on g(t) and therefore represents the excess work 
[29, 33, 34]. Since 'I'o(—t) = 'I>o(t) [see Eq. (8)], we obtain 


lAexc =/ du [ du' g{u)'i>o{T{u-u'))g{u'). 


( 12 ) 
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where g{u) and g{u') denote the derivatives with respect 
to u= {t — tQ)/T and u' = {t' — tQ)/T, respectively. 

In summary, linear response expresses the total work 
as a sum of two contributions. One is independent of 
the particular process and is identical to what would 
be obtained in the quasistatic limit, i.e., the quasistatic 
work. The other is therefore interpreted as the additional 
amount of energy that the external agent has to pump 
into the system in a finite-time process. What we call ex¬ 
cess work here is therefore defined as Wg^c = W — Wg ^, 
where Wgs is the quasistatic work. Thus, we expect that 
Eq. (12) goes to zero asymptotically as the quasistatic 
limit is approached. 


III. EXCESS WORK FOR AN EXACTLY 
SOLVABLE MODEL 


We shall apply now the expression of Eq. (12) to 
the one-dimensional harmonic oscillator, which is a com¬ 
pletely solvable model that allows us to check the accu¬ 
racy of the linear response expression of the excess work. 
Therefore, we will consider that the dynamics of our sys¬ 
tem of interest is given by the following time-dependent 
Hamiltonian: 

H[A(t)] = ^ + A(t)y, (13) 

which can model, for example, the motion of a colloidal 
particle in an optical trap [38]. 

From the solution of Hamilton’s equations for A = Aq, 
we can calculate the relaxation function exactly. Accord¬ 
ing to Eq. (8), we first need to obtain the correlation 
function C'o(t), 


a:^(0) a:^(t) \ 2 cos^ (wot) + 1 


where Wq = Aq since the mass was set equal to one in Eq. 
(13). The next step is to calculate the constant C. From 
Eqs. (9) and (14), we obtain C = 2/(4/3^Aq). Finally, the 
relaxation function reads 


'l^o(t) = = 'I^o(O) cos(2..ot), (15) 

where 'I'o(O) = (,5/2)((a:^(0)/4)o—(a:^(0)/2)Q). Although 
it is a bit misleading to call Eq. (15) a relaxation function, 
we will see next that it leads to a reasonable thermody¬ 
namic behavior of Wgxc [33, 34]. 

Substituting Eq. (15) into (12) and using the linear 
protocol g{t) = {t — to)/T, we obtain 



FIG. 2. (Color online) Comparison between numerical cal¬ 
culations (blue circles) and Eq. (16) (dashed line) for (a) 
JA/Ao = 0.1 and (b) JA/Aq = 0.5. We used 10® initial con¬ 
ditions to calculate the work for each value of the switching 
time r. 


simulations. The agreement is very good for 5A/Ao = 
0.1. Nevertheless, our linear response expression already 
deviates considerably for AA/Aq = 0.5. 

A very striking prediction of Eq. (16) is that Wexc can 
be zero for specific finite values of r. In other words, 
there are finite switching times for which the total work 
is equal to the quasistatic value. These particular values 
of r can be obtained directly from Eq. (16): whenever 
coqt = Ztt, with I integer, we have Wgxc = 0. This means 
that already for r equal to half of the natural period of 
oscillations, 2ttIujq, the system can be driven as if the pro¬ 
cess was a quasistatic one. The possibility of achieving 
the quasistatic value of the work performed in a finite¬ 
time process was pointed out before in Ref. [32], though 
without addressing the dependence of the excess work on 
the switching time. In Fig. 2, the numerical value of Wexc 
is obtained after subtracting from W the exact value of 
the quasistatic work Wgs , 


Wexcir) 


1 /<5AA ^ sin^(a;or) 

8/3 VAo j (wot)2 


(16) 



(17) 


This expression goes to zero in the quasistatic limit, r —>■ 
oo, and has its maximum value when t —> 0. Figure 
2 shows a comparison between Eq. (16) and numerical 


More details about the previous expression for Wgs can 
be found in Appendix B. Finally, we point out that the 
inset of Fig. 2(b) indicates that the minima of Wgxc are 
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FIG. 3. (Color online) Comparison between numerical cal¬ 
culations (dashed line) and Eq. (16) (blue circles) for the 
(a) quadratic, g{t) = {t — to/r)^, and the (b) exponential, 
p(t) = (1 — — e~^), protocols. In both cases, 

5X/Xo ~ 0.1. We used 10® initial conditions to calculate the 
work for each value of the switching time r. 


lifted as (5A/Ao increases. This effect is clearly beyond 
our linear response approach. 


These results tell us that a simple linear protocol can 
be optimal if we choose the value of r carefully. Be¬ 
sides, our linear response expression predicts that this 
happens only for specific finite values of t. Among the 
many interesting questions that arise from these remarks, 
we will focus now on the following: Is this a special 
feature of the linear driving? To investigate that, we 
will compare analytical and numerical results for differ¬ 
ent nonlinear protocols. Let us consider, for instance, 
a quadratic, g{t) = ((t — and an exponential, 

g{t) = (1 — e“(‘“*“^/'^)/(l — e“^), protocol. The results 
are shown in Fig. 3. The agreement is again very good 
and Wexc goes to zero in the quasistatic limit. In contrast 
to what happens for the linear driving, our analytical re¬ 
sults predict that Wexc never vanishes in finite time for 
the nonlinear protocols considered. This can be explicitly 
checked, for instance, for the quadratic protocol, whose 



FIG. 4. Examples of the family of protocols given by Eq. (19) 
for K = 0 (solid line), k = 4 (dashed line), and k = 2 (dotted 
line). We fixed a equal to one. 


expression for Wexc reads 


w 1 


8/3 \Ao 

[(wor)2 - woT sin (2a;oT) -|- sin^ (wqt)] 

"-^ 7)3 - 


Nevertheless, Fig. 3 shows that Wexc does have finite¬ 
time minima in these cases. 

It has been shown in the literature [32] that there exists 
indeed a highly degenerate family of finite-time nonlinear 
protocols for which the work performed is equal to the 
quasistatic one. In what follows, we will show how to ob¬ 
tain such protocols analytically from our linear response 
approach. 

First, we use the fact that, for g[t) = {t — to)/T, we do 
observe zeros of Wexc in finite time. Second, it can be 
easily verified that the Fourier series of the linear proto¬ 
col has no cosine coefficients in the interval to < t < tf. 
Therefore, we wonder what happens to Wexc if we per¬ 
form a protocol given by a linear part plus a sine func¬ 
tion such that the boundary conditions g{to) = 0 and 
g(tf) = 1 are preserved. In other words, we ask ourselves 
whether the protocol 

g{t) = - —+asm^KTr- -, (19) 


where k is an integer and a is an arbitrary real number, 
leads to zeros of Wexc- Some examples of these functions 
can be seen in Fig. 4. 

We show in Fig. 5 the comparison between analyti¬ 
cal and numerical calculation of Wexc using the protocol 
given by Eq. (19). These results show two additional 
features compared to those in Fig. 2: the first zero of 
Wexc occurs at shorter times and there is a constructive 
resonance for a specific value of r. These features can be 
better explained if, using Eq. (15), we rewrite Eq. (12) 
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FIG. 5. (Color online) Wej:c for the protocols (a) g(t) = (t — 
Io)/'r + sin(27r(t — to)/r) and (b) g(t) = (t — to)/r + sin(47r(t — 
to)/T). Analytical and numerical results are represented by 
dashed lines and blue circles, respectively. We used SX/Xo = 
0.1 and 10® initial conditions. 


as 

^^exc CX 

fl 

/ dug{u) cos (2ujotu) 


dug(u) sin(2wo™) 


Therefore, to have Wexc = 0, we demand that 


( 20 ) 


dug(u) cos (2u!otu) = 0 , 


dug{u) sin {2ujqtu) = 0 , 


(21a) 

(21b) 


or, equivalently, 


z = I dug{u) e 


2iujQTU _ 


= 0 . 


(21c) 


By inserting the protocol of Eq. (19) in the expression 
for z and solving the integral, we obtain 


z = 2iujoT 
+ 2ujot 


1 — cos(2wot) [1 — cos(K7r) cos(2a;oT)] 

+ aiTK- 


{2ujQTy 


(2ujot)^ - ( K 7 r )2 


sin(2u;oT) cos(K7r) sin(2wo''‘) 

—'/ -— H” CLTT _ 

(2wot) 2 (2a;or)2 - {mry 


( 22 ) 


We conclude that the real and imaginary parts of z can 
be zero simultaneously if n is even and 2ujqt = 27r/, with 
I integer, independently of the value of a. This condition 
predicts that the first zero would occur for ujqt = tt unless 
K = 2. However, we see in Fig. 5 that the position of the 
very first minima of Wexc does not follow this prediction. 
This is so because there is a second kind of zero that does 
depend on the value of a. As before, we demand that real 
and imaginary parts of z vanish, but now for the same 
value of a. We obtain then 


LUqT 


(K7r/2) 

(1 + KTTO)^/^ 


(23) 


For K = 2 and a = 1, Eq. (23) leads to loqt « 1.2, in 
agreement with Fig. 5(a). For k = A and a = 1, there 
are two zeros before the resonant peak, with the first one 
at ujqt Ri 1.7, due to the value of a, and the second one 
at ujqt = TT. Hence, we can obtain zeros of Wexc at arbi¬ 
trarily short times by choosing the values of a appropri¬ 
ately. Nevertheless, for a negative, we are limited by the 
square root in Eq. (23). Figure 6(b) shows what happens 
to Wexc when we perform the protocol of Eq. (19) with 
a = — 1 and k = 2. As opposed to Fig. 5(a), the zeros 
of Wexc in Fig- 6(b) do not depend on a. Although the 
results in Fig. 6(a) show very pronounced minima, linear 
response predicts that there are no finite-time zeros of 
Wexc in this case. 

This analysis of Wexc has interesting consequences. If 
we add to the protocol presented in Eq. (19) an arbitrary 
number of sine functions with arbitrary coefficients and 
even values of k, we still obtain zeros when wqt = hr. 
This sum of sinusoidal terms can be understood as the 
Fourier series of a function whose values at t = to nnd 
t = tf are zero and which is odd with respect to a change 
of {t — to) by T — {t — to). This property is illustrated 
in Fig. 7. Therefore, we conclude that any function that 
vanishes at to and tf and is odd with respect to time 
reversal leads to the above-mentioned zeros of Wexc when 
added to the linear protocol. For instance, the following 
family of polynomials: 


fk{t) = € 


1 - 2((t - to)/T) 

22fc+l 


t - to 


\ 2/c+1 

2/ 


(24) 

where k is an integer and e = ±1, has such property, as 
illustrated in Fig. 8. For k = 1, the sinusoidal term of 
Eq. (19) is not odd with respect to time reversal. 


IV. NONLINEAR SYSTEMS 

The linear response expression (12) depends strongly 
on the behavior of a specific autocorrelation function. 
For the system (13), this function oscillates indefinitely 
when the dynamics is time independent. This is a special 
feature of linear systems whose observables have frequen¬ 
cies of oscillations which are independent of the energy. 
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FIG. 6. (Color online) Wexc for the protocol (19) with (a) 
K = 1 and a = 1 and (b) k = 2 and a = — 1. Analytical and 
numerical results are represented by dashed lines and blue 
circles, respectively. We used S\/Xo = 0.1 and 10® initial 
conditions. 


In contrast, the dynamics of nonlinear systems, as 

’^['^(0] = ^ +1 (25a) 

and 

^[AW] = y + A(i)y, (25b) 

present frequencies of oscillations which are energy de¬ 
pendent. Thus, correlation functions are oscillatory only 
when initial conditions are sampled from a single energy 
shell. Otherwise, a decay is observed due to the incom¬ 
mensurability of the superposed oscillations from differ¬ 
ent energy shells. Figure (9) shows an example of this 
for the system given in Eq. (25a). 

In this section, we argue that the previous analysis 
of Wexc can be extended to one-dimensional anharmonic 
oscillators. In particular, we want to investigate the be¬ 
havior of Wexc when systems (25a) and (25b) are driven 
by the protocols discussed in the previous section. As 
shown in Fig. 9, the relevant correlation function of sys¬ 
tem (25a) has well-defined oscillations for short times. 
For the sake of clarity, let us assume for a moment that 
these oscillations last indefinitely with a period Ta- In 



FIG. 7. (Golor online) Sinusoidal terms, sin (K7r(t — to)/r), 
for K = 2 (solid line) and k = 4 (dotted line), and 
— sin («:7r(t — to)/T) for k, — 4 (dashed line). 



FIG. 8. Examples of the polynomials given by Eq. (24) for 
k = 1 and e = -|-1 (dotted line), k — 2 and e = -|-1 (dashed 
line), fc = 1 and e = —1 (solid line), and k = 2 and e = —1 
(dash-dotted line). 


this case, we can replace the actual d'o(t) by its Fourier 
series in the interval [—Ta/ 2, T/i/2], 


OO 

= d<o(0) an cos{nujAt) , (26) 

71 — 1 


recalling that = lE'o(i)- The a„ are the Fourier 

coefficients and loa = 2-k/Ta [according to Eq. (8), d'o(t) 
would oscillate around zero implying that the coefficient 
oq is zero]. Substituting the previous expression into 
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FIG. 9. Numerical calculation of (2;^(0)2:^(t))o/4 for the os¬ 
cillator (25a) with time-independent A equal to Aq. Initial 
conditions were sampled according to a Boltzmann-Gibbs dis¬ 
tribution with Ao = 1. The period of oscillations Ta is ap¬ 
proximately equal to 1.74 in arbitrary units. We used 10® 
initial conditions. 


Eq. (12), we have 
= ^'I'O(O) 

CJO 

X flfc du du' g{u) cos{nLo at{ u — u'))g{u'). 

n=i Jo Jo 


(27) 


Analogously to Sec. Ill, Eq. (27) can be written as 




du g(u) cos{nuJATu) 


-\- du g{u)sin{nu!ATu) 


(28) 


where A„ = ((i5A)^/2)'I'o(0) Un- This expression for the 
excess work vanishes if, for instance, each term of the sum 
is zero for the same value of r. To verify this possibility, 
we check under what conditions 


= 



du g{u) 


^inuiATU _ Q 


(29) 


Using the protocol of Eq. (19), the quantity za reads 


2A = 


inoJAT 


1 — cosiriuJAT) 


nOJAT 


{jiloatY 

sin(nwA'r) 


{uujatY 


aiTK 


(1 — cos(K7r) cos{nu}AT)) 
{nujATY ~ (K7r)2 
cos(K7r) sin(nwA'r) 


{tvjjatY ~ [k-tY 


(30) 


Thus, for K even, there are zeros of za whenever hloat = 
27rZ„, with In an integer, except for 2Z„ = k. Therefore, 
for K = 0, the smallest value of ujat that provides a zero 


of ZA for all modes simultaneously is 27r. When k = 2, 
this zero is forbidden by the denominators in Eq. (30) 
and ujat = tt is the first zero. 

This prediction of the first minimum of Wexc is in very 
good agreement with the numerical calculations shown 
in Fig. 10. However, it is based on a wrong assump¬ 
tion about the behavior of 4'o(i)- Introducing a small 
damping of oscillations, the Fourier transform of d'o(t) 
changes from a delta-like peak at uja to a peak with a 
small width whose position is very close (but not exactly 
equal) to uja- Therefore, instead of the representation 
given by Eq. (26), we would have the following one: 

^'o(i) = ^o(0)y^y dwx(w)cos(wt), (31) 

where x(w) is the cosine Fourier transform of 
'I'o(t)/'I'o(0). It is not hard to see that by plugging the 
expression (31) into Eq. (12), we obtain an expression 
similar to (28) with the sum replaced by the integral 
over w and the coefficient A„ replaced by (y/2/7r)x(a;). 
In this case, we have a continuous (but small) interval 
of frequencies that contribute to W^xc instead of the dis¬ 
crete values nujA- Then, we can think of a z^(w) given 
exactly by Eq. (30) with uuja replaced by w. To observe 
a zero in the excess work, za(uj) would have to be zero for 
all a; in a small vicinity of uja = uja + i 5 (<5 represents the 
small shift of the peak due to damping). However, due to 
the incommensurability of such frequencies, if za{uja) is 
exactly zero, it is certainly nonzero for w around uja- In 
other words, if za{uja) = 0, then z^(a;) ~ 0 for w ~ Qa, 
and a minimum of Wexc arises. Therefore, our linear re¬ 
sponse approach predicts that Wexc does not vanish in 
finite time for the family of protocols (19) when 'I'o(0 
has damped oscillations. 

As mentioned before. Fig. 10 illustrates the preceding 
discussion. The value of Wqs was obtained analytically in 
Appendix B. Figure 10(a) shows a numerical calculation 
of Wexc for the system given by expression (25a) using 
the linear protocol. We kept 5X/Xo = 0.1 so that our 
approach based on linear response theory is still valid. 
This numerical result shows that Wexc has indeed minima 
and the position of the first one is very close to ujat = 

The numerical calculation of Wexc for the protocol given 
by Eq. (19) is shown in Fig. 10(b). It shows the same 
features we observe in Fig. 5 for the harmonic oscillator, 
including a minimum at very short time scales for ujat ~ 
1.2. As discussed in Sec. HI, this minimum is related to 
the value of a. For the system (25b), the numerical result 
is shown in Fig. 11 using the linear protocol. Although 
the excess work also approaches zero as for a finite r, we 
cannot distinguish between a minimum and a monotonic 
decay. 

This analysis of Wexc for the system (25a) shows how 
determinant is the behavior of the relaxation function. If 
correlations decay sufficiently fast, a crossover is observed 
in the behavior of Wexc- This can be easily verified using 
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FIG. 10. (Color online) Numerical calculation of W^^c. for the 
anharmonic oscillator (25a) using (a) a linear protocol and 
(b) the protocol (19) with k = 2 and a = 1. In both cases, 
5A/Ao = 0.1 and we used 10® initial conditions. The vertical 
dotted lines indicate the analytical prediction, (a) ujat = 27r 
and (b) ujat ~ 1.2, of the first minimum. 


the phenomenological expression 



FIG. 11. (Color online) Numerical calculation of We^c for 
the anharmonic oscillator (25b) using the linear protocol and 
(5A/Ao = 0.1. We used 10® initial conditions. The inset shows 
a small region around r = 5.2. 
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4'(t) = T'o(O) e “1*1 \ cos{uJAt) + ^ sin(a;A|t|)j , (32) 

where a and uja denote the decay rate and the frequency 
of oscillations, respectively. The excess work obtained for 
the linear protocol using Eqs. (32) and (12) is shown in 
Fig. 12 for different ratios of ajujA- As the decay rate in¬ 
creases, the minima disappear and Wexc varies monoton- 
ically with r. We have verified the same sort of crossover 
when correlations decay as a power law. 

V. EXCESS WORK AND THE ADIABATIC 
INVARIANT 

Motivated by the results of the previous sections, we 
argue here how the excess work is connected to an impor¬ 
tant quantity of time-dependent Hamiltonian systems, 
namely, the adiabatic invariant. Adiabatic invariants are 
approximate constants of motion of time-dependent sys¬ 
tems perturbed by slowly varying parameters [39, 40]. 
Generally, the perfect conservation of an adiabatic invari¬ 


FIG. 12. (Color online) Excess work, Wexc = 
2II4a;c/(((5A)^'ko(0)), obtained from expression (32) for 
ajuJA = 0.01 (red dotted line), 0.1 (blue dashed line), and 
0.3 (black solid line) using the linear protocol. 


ant fl is reached only in the quasistatic limit. Besides, 
for systems with one degree of freedom, it is possible to 
estimate analytically the conservation of H with high ac¬ 
curacy [41]. For this class of systems, the results obtained 
previously show that it is indeed possible to have W ap¬ 
proximately or even exactly equal to Wqs for finite values 
of switching time. Thus, it is reasonable to expect that 
the H is conserved (or almost conserved) at exactly those 
values of r for which Wexc vanishes (or has minima). 

In what follows, we investigate the behavior of AH = 
H(t/) — H(to) for systems given by expressions (13) and 
(25a). The adiabatic invariant H of both oscillators is 
equal to the area enclosed by the energy shell and, ac¬ 
cording to Appendix A, it can be expressed in terms of 
the energy and A only. Therefore, for fixed values of Aq 
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FIG. 13. (Color online) Numerical calculations of VFea;c (blue 
dots) and (dashed line) for the harmonic oscillator (13) 
using the linear protocol. We fixed (5A/Ao = 0.1 and used 10® 
initial conditions. 

and A/, AH depends exclusively on the initial and final 
energies as we change r. It is shown in Ref. [41] that if we 
sample initial conditions with the same value of energy, 
say Eq, and evolve the equations of motion of system 
(13) for an arbitrary time interval r using a given proto¬ 
col A(t), the distribution of final energies Ei is such that 


{Ei-Ei)^’^+^ =0, (33a) 

(E, - ^7)2- = (^{E,-W,ry , (33b) 

where m is an integer and the overbar denotes an average 
on the distribution of Ei. Equations (33) show that any 
moment of this distribution can be written in terms of the 
averaged energy Ei and the variance fj? = {Ei — Ei)^. 
Hence, there is an exact conservation of the adiabatic in¬ 
variant whenever = 0 because the system would have 
evolved from a single energy shell to another. In the re¬ 
mainder of this section, we will compare numerical results 
of and for both systems mentioned previously. 

We see in Fig. 13 a clear agreement between the be¬ 
haviors of Wexc and jj? for the harmonic oscillator. This 
result was obtained using a linear protocol. Furthermore, 
it is shown in Fig. 14 that agreement is also very good for 
protocols of the family (19). These results suggest that 
the finite-time zeros of Wexc imply the conservation of 
H. Conversely, we would like to have a proof that every 
time H is conserved in finite time, W^xc vanishes. For 
the moment we only have numerical evidence that when 
the Hamiltonian (13) is driven by the family of protocols 
(19), H is conserved whenever Wexc vanishes in finite 
time. We have observed this no matter the initial energy 
shell we start. 

The connection between the adiabatic invariant and 
thermodynamic work was not mentioned in previous 
works about optimal paths in thermally isolated systems. 



LJoT 

FIG. 14. (Golor online) Gomparison between the analytical 
prediction of Wexc (dashed line) given by (12) and the numer¬ 
ical calculation of for the harmonic oscillator (13) using 
the protocol g(t) = (t — to)/T -|- sin(27r(t — to)/T). We fixed 
SX/Xo ~ 0.1 and used 10® initial conditions. 



FIG. 15. (Golor online) Comparison between numerical cal¬ 
culations of Wexc (squares) and (dashed line) for the an- 
harmonic oscillator (25a) using the linear protocol. We fixed 
(5A/Ao = 0.1 and used 10® initial conditions. 


Figure 15 suggests that such relation also exists for the 
anharmonic potential (25a). However, the zeros of 
do not imply the conservation of H in this case because 
Eqs. (33) do not apply [42, 43], i.e., the vanishing of 
does not imply that all higher-order moments vanish 
too. This can be verified numerically, constructing the 
distribution of Ei for the values of r where vanishes. 
Figure 16 shows that this distribution does not have a 
single peak, indicating that at the end of the protocol 
there is not just one single value of H. 
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El 


FIG. 16. (Color online) Energy distribution at the end of the 
linear protocol for system (25a). We have chosen ljJat such 
that jj,^ has its first zero at this value (see Fig. 15). Initial 
conditions were sampled from a single energy shell. 


VI. DISCUSSION 

As mentioned before, the existence of finite-time pro¬ 
cesses leading to Wexc = 0 was first reported in Ref. [32] 
for thermally isolated systems. There, the authors claim 
the existence of highly degenerate protocols for which 
the work performed is equal to the quasistatic value. In 
the regime described by our linear response approach, we 
were able not only to confirm the existence of such pro¬ 
tocols but to show that they must obey a certain sym¬ 
metry. According to Sec. Ill, the optimal finite-time 
protocols we found are composed of two parts: a lin¬ 
ear protocol plus a function which is odd with respect to 
time reversal. This symmetry embraces a much larger 
class of protocols than those of Ref. [32]. On the other 
hand, we have shown that these optimal protocols never 
lead to Wexc = 0 for one-dimensional anharmonic oscilla¬ 
tors. Instead, Wexc has minima when such protocols are 
performed. The problem of Wexc = 0 in finite time for 
nonlinear systems was addressed in Ref. [32] and has re¬ 
mained inconclusive. Unfortunately, our contribution to 
this problem is very restricted: we have shown that the 
family (19) does not lead to finite-time zeros of Wexc- 
Thus, our results do not exclude the existence of other 
families of optimal protocols. 

The fact that Wexc, as a function of r, can have finite¬ 
time zeros or minima contradicts the usual intuition of a 
monotonic decay as the quasistatic limit is approached. 
At a first glance, this may be wrongly taken as an ex¬ 
clusive feature of systems with few degrees of freedom. 
However, Eq. (12) tells us that what really matters is the 
behavior of the relaxation function. Therefore, we can in¬ 
fer from the analysis of the simple models presented here 
that this apparently peculiar thermodynamic behavior 
might also show up in the thermodynamic limit. This is 
the great advantage of our phenomenological approach. 


Indeed, it is well known that Eq. (32) can describe the 
decay of correlations of a large class of macroscopic sys¬ 
tems [37] . Although we leave for a future work the study 
of more complex systems, we would like to briefly out¬ 
line how the physics of Eq. (12) provides approximate 
solutions of the optimization problem in this case. If the 
relaxation function is known from computer simulations, 
then Eq. (12) tells us that we must find a g{u) such that 
the surface generated by II'o('7‘(m — u'))g{u)g{u') has the 
minimal volume inside the integration domain. In other 
words, from the knowledge of 'l'o(t), it is possible to find 
approximate optimal solutions by geometric inspection. 

Another aspect of the nonmonotonic decay of Wexc is 
worth mentioning. As the switching time goes to zero, 
our results approach the value PiWjp — Wqs) no matter 
the protocol we use. The quantity Wjp represents the 
work performed when A is suddenly switched from Aq 
to A/ and its value is simply (U(A/) — U(Ao))o, where 
V (A) is the potential energy. Thereby, this is the fastest 
protocol we can perform. One would expect then that 
Wexc is maximum when r —>■ 0. However, this is not 
what we have observed. According to Figs. 5, 6, and 10, 
there are indeed finite-time peaks of Wexc whose values 
are much larger than Wexc{T —t 0). 

We also want to point out that definitions (21c) and 
(30) are essentially the rapidity parameter appearing in 
the study of quantum work distributions of the thermally 
isolated harmonic oscillator [32, 44, 45]. From what was 
shown here, we believe that za probably plays an impor¬ 
tant role in the statistics of quantum work of anharmonic 
oscillators [46]. 

As a last remark, we want to mention the relation be¬ 
tween the excess work and the conservation of the adi¬ 
abatic invariant. Our results suggest that, in general, 
finite-time zeros (or minima) of Wexc imply the conser¬ 
vation (or almost conservation) of il. This relation de¬ 
serves a more careful analysis due to its potential use¬ 
fulness in the search for optimal paths. Besides, this 
has interesting implications to the adiabatic switching 
method proposed by Watanabe and Reinhardt [5] to es¬ 
timate entropy and free-energy differences: whenever H 
is almost conserved in finite time, then this method will 
also provide good estimates in finite time. As discussed 
by the authors in Ref. [5], the problem of course is how 
to find, in general, the switching protocol that does the 
job. Our approach suggests that although the existence 
of such optimal paths depends very much on the dynam¬ 
ics of the system, it could be inferred from the behavior 
of the corresponding relaxation function (see Fig. 12). 


VII. CONCLUSIONS 

In summary, we have shown that within linear re¬ 
sponse, there are highly degenerate protocols which lead 
to finite-time zeros or minima of the so-called excess 
work, Wexc, on a thermally isolated system. This quan¬ 
tity was defined as the amount of energy the external 
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agent has to pump into the system in addition to the 
quasistatic work. Therefore, every time the excess work 
vanishes in finite time, the total work is equal to the qua¬ 
sistatic value. According to our approach, the family of 
optimal protocols must be composed of a linear part plus 
a function which is odd with respect to time reversal. 

Our analytical and numerical results have shown a 
counterintuitive behavior of the excess work as a func¬ 
tion of the switching time, namely, a nonmonotonic de¬ 
cay as the process becomes slower. Although obtained 
for small systems, we claim that this behavior exists in 
macroscopic systems as well. Our argument relies on the 
expression for the excess work based on the relaxation 
function. We have shown that, for weak enough decay 
of correlations, this effect must be present no matter the 
size of the system. In other words, the only requirement 
is that the driving force dH / d\ has a sufficiently oscilla¬ 
tory autocorrelation function. 

Finally, the relation between finite-time zeros (or min¬ 
ima) of Wexc and the conservation of the adiabatic invari¬ 
ant n suggests that there may exist an interesting and 
useful connection between optimal finite-time processes 
and shortcuts to adiabaticity [47, 48]. 
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1. Harmonic Oscillator 


For the Hamiltonian (13), the transformation men¬ 
tioned before reads 


p = cos (0), X = \/2'H/\ sin {8) 
and its Jacobian is given by 


V2'H cos(0) s\n{8)/\/2'H 
— shi{9)J ^ cos{6)/y/2'H\ 




Hence, the adiabatic invariant reads 

2nE 


1 


n{E,x) = 


Vx ■ 


(A2) 


= ^- (A3) 


(A4) 


2. Anharmonic Oscillator I 


Considering now the Hamiltonian (25a), the phase- 
space parametrization is given by 

= f ' cos^A (0), (A5) 


p = ypXH. sin (0), a; = 


A / 


where 0 < 0 < 7r/2 and its Jacobian reads 


J{8,U) = 


y/mcos{e) suY{e)/y/m 


1/2^ M/4 sin(e) H. 
A > cosi/2(e) 

^-1/4 




(S) 


25/4A1/4 cos "/^(6»)] 


27/4 Xl/i 


Performing the integrals, we obtain 


(A6) 


(A7) 


where K im) is the incomplete elliptic integral of the first 
kind given by 

f.7r/2 

(A8) 


K{m) = [ 
Jo 


d(j) 


0 — m sin^ (0) 


Appendix A: Adiabatic invariant 


3. Anharmonic Oscillator II 


In this Appendix, we obtain the adiabatic invariant 
H(i?, A) for the systems considered previously. Since it 
is the area enclosed by the energy shell [40], H(A, A) can 
be calculated as follows: 


H(i?, A) = J dx dpQ{E —'H{x,p] X)) 

pE p2-K 

= dn dOJiO.'H), (Al) 
Jo Jo 


where is the Jacobian of the transformation {x,p) —/ 

{o.n)- 


The phase-space parametrization for Hamiltonian 
(25b) is given by 

1/6 


p = y/^Hsin (8 ), 




cos3/3 (0), (A9) 


where 0 < 8 < 12. After calculating the Jacobian 

sin(0)/-\/^ 


J{8,'H) = 


\/‘2'H cos{d) 

21/6-H3/3 sin(e) W-3/6 l/3(n\ 

35/6Ai/a cos2/3(e) 63/3Ai/6 ) 

-W-l/S r 

sin^(0) cos“3/3(0) + cos‘^^^{9) 


21/335/6^^6 L’ 


(AlO) 
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we obtain 


22/^3V6^r(7/6) 2/3 

AV6 r(2/3)^ ’ 

where r(a;) is the gamma function defined as 


which implies 


(All) 





3. Anharmonic Oscillator II 


(B5) 


r(x) 





(A12) 


Finally, for the anharmonic oscillator (25b), the conser¬ 
vation of its adiabatic invariant (All) along a quasistatic 
process provides 


Appendix B: Calculation of quasistatic work 

In this Appendix, we derive exact analytical expres¬ 
sions for the quasistatic work. This is the quantity we 
have to subtract from the numerical value of the mean 
work W to obtain the excess work Wexc- Since the system 
is thermally isolated during its time evolution, we have 
from the first law of thermodynamics that W = AU, 
where U is the internal energy. Hence, in the quasistatic 
limit, 


Wg, = AU = {Ef) - (A,). (Bl) 

where is the energy obtained via the conservation 
of the adiabatic invariant as a function of the initial 
energy, and the initial and final values of the con¬ 
trol parameter A. The brackets (.) denote an average 
on a Boltzmann-Gibbs distribution since the system was 
initially in equilibrium with a heat bath. Besides, this 
distribution is taken with A = Ai, which is the initial 
value of A. 


1. Harmonic Oscillator 

The adiabatic invariant of the harmonic oscillator (13) 
is given by Eq. (A4). Thus, for a quasistatic process, the 
conservation of H(A, A) yields 



which implies 


Wgs 


1 



1 


(B2) 


(B3) 


2. Anharmonic Oscillator I 

Considering now a quasistatic process performed on 
the anharmonic oscillator (25a), the conservation of 
fi(E, A) given by (A7) yields 


Ef = E, 



(B4) 


Ef = E, 



which implies 


Wgs 



(B6) 


(B7) 


4. Linear Response expression of Wg. 


We can now compare the linear response expression for 
the quasistatic work. 


W, 




\ 9A /f 


^o(O), 


with the exact results derived in Appendix A. 
For the harmonic oscillator (13), we have 


d'H \ lx" 


1 


2/3 Ao 


(B8) 


(B9) 


and 


4-0 = «'o(0)-x“ 

T 


4^A, 


2 ’ 
0 


(BIO) 


since Xcf = 0. 

Hence, Eq. (B8) reads 


1 


IT = - 

qs 


sx 1 rsx 

2Ao 8 I Ao 


(Bll) 


The expansion of the exact result (B3) for SX/Xq 1 
reads 


_ 1 


E] -1 




XoJ 


SX 1 /SX 


2Ao 8 \ Ao 


0(3), 


(B12) 


which is equal to (Bll) up to second order in SX/Xq. 
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